Chapter 14 Geospatial R - Raster Hydro Analyses

14.1 Introduction

The following activity is available as a template github repository at the following link:

For more:

Load necessary packages and data

#install.packages("whitebox", repos="http://R-Forge.R-project.org")

library(tidyverse)
library(raster)
library(sf)
library(whitebox)
library(tmap)
whitebox::wbt_init()

theme_set(theme_classic())

tmap_mode("plot")
# brush <- raster("McDonaldHollowDEM/brushDEMsmall.tif")
# 
# wbt_resample("McDonaldHollowDEM/brushDEMsmall.tif",
#              "McDonaldHollowDEM/brushDEMsm_5m.tif",
#              cell_size = 0.00003925595)
# 
# brush <- raster("McDonaldHollowDEM/brushDEMsm_5m.tif")
# 
# tm_shape(brush)+
#   tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
#   tm_scale_bar()
tmap_mode("view")
## tmap mode set to interactive viewing
dem <- raster("McDonaldHollowDEM/brushDEMsm_5m.tif")
 

dem[dem < 1500] <- NA

tm_shape(dem)+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape dem unknown. Long lat (epsg 4326)
## coordinates assumed.

Hillshade for visualization

wbt_hillshade("McDonaldHollowDEM/brushDEMsm_5m.tif",
              "McDonaldHollowDEM/brush_hillshade.tif")
## [1] "hillshade - Elapsed Time (excluding I/O): 0.12s"
hillshade <- raster("McDonaldHollowDEM/brush_hillshade.tif")

tm_shape(hillshade)+
  tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
  tm_scale_bar()
## Warning: Currect projection of shape hillshade unknown. Long lat (epsg 4326)
## coordinates assumed.

Prep DEM for hydro analysis

wbt_fill_single_cell_pits("McDonaldHollowDEM/brushDEMsm_5m.tif",
                     "McDonaldHollowDEM/bmstationdem_filled.tif")
## [1] "fill_single_cell_pits - Elapsed Time (excluding I/O): 0.4s"
wbt_breach_depressions_least_cost("McDonaldHollowDEM/bmstationdem_filled.tif",
                     "McDonaldHollowDEM/bmstationdem_filled_breached.tif",
                     dist = 5,
                     fill = TRUE)
## [1] "breach_depressions_least_cost - Elapsed Time (excluding I/O): 0.52s"
filled_breached <- raster("McDonaldHollowDEM/bmstationdem_filled_breached.tif")

## What did this do?
difference <- dem - filled_breached

difference[difference == 0] <- NA

tm_shape(hillshade)+
  tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
  tm_scale_bar()+
tm_shape(difference)+
  tm_raster(style = "cont",legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape hillshade unknown. Long lat (epsg 4326)
## coordinates assumed.
## Warning: Currect projection of shape difference unknown. Long lat (epsg 4326)
## coordinates assumed.
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
wbt_slope("McDonaldHollowDEM/brushDEMsm_5m.tif",
          output = "McDonaldHollowDEM/demslope.tif",
          units = "degrees")
## [1] "slope - Elapsed Time (excluding I/O): 0.5s"
slope <- raster("McDonaldHollowDEM/demslope.tif")

tm_shape(slope)+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape slope unknown. Long lat (epsg 4326)
## coordinates assumed.
wbt_d8_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
                         "McDonaldHollowDEM/D8FA.tif",
                         out_type = "specific contributing area")
## [1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 0.24s"
d8 <- raster("McDonaldHollowDEM/D8FA.tif")

tm_shape(log(d8))+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape log(d8) unknown. Long lat (epsg 4326)
## coordinates assumed.
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
wbt_d_inf_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
                         "McDonaldHollowDEM/DinfFA.tif",
                         out_type = "specific contributing area")
## [1] "d_inf_flow_accumulation - Elapsed Time (excluding I/O): 0.37s"
dinf <- raster("McDonaldHollowDEM/DinfFA.tif")

tm_shape(log(dinf))+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape log(dinf) unknown. Long lat (epsg 4326)
## coordinates assumed.
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
wbt_wetness_index(sca = "McDonaldHollowDEM/D8FA.tif",
          slope = "McDonaldHollowDEM/demslope.tif",
          output = "McDonaldHollowDEM/TWI.tif")
## [1] "wetness_index - Elapsed Time (excluding I/O): 0.4s"
twi <- raster("McDonaldHollowDEM/TWI.tif")

tm_shape(twi)+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape twi unknown. Long lat (epsg 4326)
## coordinates assumed.
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
wbt_downslope_index("McDonaldHollowDEM/brushDEMsm_5m.tif",
                      "McDonaldHollowDEM/TWId.tif")
## [1] "downslope_index - **********************************************************************************"
twid <- raster("McDonaldHollowDEM/TWId.tif")
twid[twid > 4000000 | twid <= 0] <- NA

tm_shape(log(twid))+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape log(twid) unknown. Long lat (epsg 4326)
## coordinates assumed.
wbt_d8_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
                            "McDonaldHollowDEM/D8FA.tif")
## [1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 0.20s"
d_inf <- raster("McDonaldHollowDEM/D8FA.tif")

tm_shape(log(d_inf))+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()
## Warning: Currect projection of shape log(d_inf) unknown. Long lat (epsg 4326)
## coordinates assumed.
wbt_extract_streams("McDonaldHollowDEM/D8FA.tif",
                    "McDonaldHollowDEM/raster_streams.tif",
                    threshold = 6000)
## [1] "extract_streams - Elapsed Time (excluding I/O): 0.4s"
wbt_d8_pointer("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
               "McDonaldHollowDEM/D8pointer.tif")
## [1] "d8_pointer - Elapsed Time (excluding I/O): 0.5s"
wbt_raster_streams_to_vector("McDonaldHollowDEM/raster_streams.tif",
                              "McDonaldHollowDEM/D8pointer.tif",
                              "McDonaldHollowDEM/streams.shp")
## [1] "raster_streams_to_vector - Elapsed Time (excluding I/O): 0.4s"
streams <- st_read("McDonaldHollowDEM/streams.shp")
## Reading layer `streams' from data source `/Volumes/GoogleDrive/My Drive/CLASSES/Hydroinformatics/Hydroinformatics_Bookdown/McDonaldHollowDEM/streams.shp' using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 2 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -80.49733 ymin: 37.23597 xmax: -80.46388 ymax: 37.25634
## CRS:            NA
tm_shape(streams)+
  tm_lines()+
  tm_scale_bar()
## Warning: Currect projection of shape streams unknown. Long-lat (WGS84) is
## assumed.
points <- read_csv("McDonaldHollowDEM/points.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   lon = col_double(),
##   lat = col_double()
## )
pointsSP <- SpatialPoints(points, proj4string = CRS("+proj=longlat +datum=WGS84"))

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(log(twid))+
  tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
  tm_scale_bar()+
tm_shape(pointsSP)+
  tm_dots(col = "red")
## Warning: Currect projection of shape log(twid) unknown. Long lat (epsg 4326)
## coordinates assumed.

twidvals <-  extract(x = twid, #raster 
                 y = pointsSP, #points
                 method = "simple")
points$twid <- twidvals

points
## # A tibble: 3 x 3
##     lon   lat    twid
##   <dbl> <dbl>   <dbl>
## 1 -80.5  37.2  89693.
## 2 -80.5  37.2 152567.
## 3 -80.5  37.2  54294.
raster_stack <- stack(twi, twid, slope, dem)

raster_values <-extract(x = raster_stack, #raster 
                 y = pointsSP, #points
                 method = "simple")

points <- cbind(points, raster_values)

points
##         lon      lat      twid       TWI      TWId demslope brushDEMsm_5m
## 1 -80.48355 37.24087  89693.31 -8.086087  89693.31 37.42480      2130.564
## 2 -80.48705 37.24570 152567.28 -9.103719 152567.28 54.68245      2301.624
## 3 -80.48477 37.24697  54293.62 -6.767529  54293.62 30.11527      2430.992